notebook.community

Edit and run 





In [2]:



    


import numpy as np
import matplotlib.pyplot as pl
import seaborn as sn
import pymilne
import scipy.integrate
%matplotlib inline



    











In [2]:



    


stI=np.load('../data/stokesI_sunspot.npy')
stV=np.load('../data/stokesV_sunspot.npy')



    











In [3]:



    


stI.shape



    


















    
Out[3]:








(512, 512, 112)

Index

Weak-field approximation for Stokes V

Strong-field approximation for Stokes V

Longitudinal magnetograph

Center-of-gravity

Unresolved fields - inclination angle

Unresolved fields - line ratio

Weak-field approximation Stokes V

Let us use the linear relation between the circular polarization profile and the derivative of the intensity profile. This amounts to solving a linear fit problem, that can be analytically obtained:

$$ V = \alpha \frac{dI}{d\lambda} $$

This is equivalent to doing the linear fit $y=ax$, which can be trivially solved by writing the $\chi^2$ merit function, differentiating it and equating to zero:

$$ \chi^2 = \sum_i (y_i-ax_i)^2 $$

We solve

$$ \frac{d \chi^2}{da} = -2\sum_i (y_i-a x_i)x_i=0, $$

with the solution

$$ a = \frac{\sum_i y_i x_i}{\sum_i x_i^2} $$ 




In [42]:



    


lambda0 = 6301.5
JUp = 1.0
JLow = 1.0
gUp = 2.5
gLow = 0.0
lambdaStart = 6300.8
lambdaStep = 0.01
nLambda = 150
wavelength = lambdaStart + np.arange(nLambda) * lambdaStep

lineInfo = np.asarray([lambda0, JUp, JLow, gUp, gLow, lambdaStart, lambdaStep])

s = pymilne.milne(nLambda, lineInfo)
BField = 100.0
BTheta = 0.0
BChi = 0.0
VMac = 0.0
damping = 0.04
B0 = 1.0
B1 = 5.0
mu = 1.0
VDop = 0.080
kl = 6.75

modelSingle = np.asarray([BField, BTheta, BChi, VMac, damping, B0, B1, VDop, kl])
stokes = s.synth(modelSingle,mu)
cont = B0+B1



    











In [43]:



    


f, ax = pl.subplots(ncols=2, nrows=1, figsize=(12,6))
ax[0].plot(wavelength, stokes[0,:] / cont)
ax[1].plot(wavelength, stokes[3,:] / cont)
ax[0].set_xlabel('Wavelength [$\AA$]')
ax[1].set_xlabel('Wavelength [$\AA$]')
ax[0].set_ylabel('I/I$_c$')
ax[1].set_ylabel('V/I$_c$')
ax[0].ticklabel_format(useOffset=False)
ax[1].ticklabel_format(useOffset=False)
pl.tight_layout()



    


















    
























In [47]:



    


geff = 0.5*(gUp+gLow) + 0.25*(gUp-gLow)*(JUp*(JUp+1)-JLow*(JLow+1))
diffI = np.gradient(stokes[0,:], lambdaStep)
alpha = -4.6686e-13*6301.5**2*geff
print "Estimated field = {0} G".format(np.sum(stokes[3,:] * diffI) / np.sum(diffI**2) / alpha)



    


















    






Estimated field = 100.520262963 G

Strong-field approximation Stokes V

In the strong field approximation, the separation of the components is just proportional to $\Delta \lambda_B$. So we only need to calculate this separation and get the field.





In [49]:



    


BField = 10500.0
BTheta = 0.0
BChi = 0.0
VMac = 0.0
damping = 0.04
B0 = 1.0
B1 = 5.0
mu = 1.0
VDop = 0.080
kl = 6.75

modelSingle = np.asarray([BField, BTheta, BChi, VMac, damping, B0, B1, VDop, kl])
stokes = s.synth(modelSingle,mu)
cont = B0+B1



    











In [50]:



    


f, ax = pl.subplots(ncols=2, nrows=1, figsize=(12,6))
ax[0].plot(wavelength, stokes[0,:] / cont)
ax[1].plot(wavelength, stokes[3,:] / cont)
ax[0].set_xlabel('Wavelength [$\AA$]')
ax[1].set_xlabel('Wavelength [$\AA$]')
ax[0].set_ylabel('I/I$_c$')
ax[1].set_ylabel('V/I$_c$')
ax[0].ticklabel_format(useOffset=False)
ax[1].ticklabel_format(useOffset=False)
pl.tight_layout()



    


















    
























In [51]:



    


maxLoc = np.argmax(stokes[3,:])
deltaL = lambda0 - wavelength[maxLoc]
print "Estimated field = {0} G".format(deltaL / (4.6686e-13*6301.**2*gUp))



    


















    






Estimated field = 10574.2712437 G

Longitudinal magnetograph





In [2]:



    


lambda0 = 6302.5
JUp = 1.0
JLow = 1.0
gUp = 2.5
gLow = 0.0
lambdaStart = 6301.8
lambdaStep = 0.03
nLambda = 50
wavelength = lambdaStart + np.arange(nLambda) * lambdaStep

lineInfo = np.asarray([lambda0, JUp, JLow, gUp, gLow, lambdaStart, lambdaStep])

s = pymilne.milne(nLambda, lineInfo)

BField = 100.0
BTheta = 0.0
BChi = 0.0
VMac = 0.0
damping = 0.04
B0 = 1.0
B1 = 5.0
mu = 1.0
VDop = 0.040
kl = 6.75

pV = np.exp(-(wavelength-lambda0-60e-3)**2 / 75e-3**2)
pQ = np.exp(-(wavelength-lambda0-120e-3)**2 / 75e-3**2)

nB = 100
B = np.linspace(0.0,3500.0,nB)

Bpar = np.zeros((2,4,nB))
Bperp = np.zeros((4,nB))
SV = np.zeros((2,4,nB))
SQ = np.zeros((4,nB))

pal = sn.color_palette()

thetas = [0, 30, 60, 85]

for j in range(4):
	for i in range(nB):
		modelSingle = np.asarray([B[i], thetas[j], 0.0, VMac, damping, B0, B1, VDop, kl])
		stokes = s.synth(modelSingle,mu)
		Bpar[0,j,i] = B[i] * np.cos(thetas[j]*np.pi/180.0)
		Bperp[j,i] = B[i] * np.sin(thetas[j]*np.pi/180.0)
		SV[0,j,i] = scipy.integrate.simps(stokes[3,:] * pV, x=wavelength) / scipy.integrate.simps(stokes[0,:] * pV, x=wavelength)
		SQ[j,i] = -scipy.integrate.simps(stokes[1,:] * pQ, x=wavelength) / scipy.integrate.simps(stokes[0,:] * pQ, x=wavelength)
		
		modelSingle = np.asarray([B[i], 180.0-thetas[j], 0.0, VMac, damping, B0, B1, VDop, kl])
		stokes = s.synth(modelSingle,mu)
		Bpar[1,j,i] = B[i] * np.cos((180.0-thetas[j])*np.pi/180.0)
		SV[1,j,i] = scipy.integrate.simps(stokes[3,:] * pV, x=wavelength) / scipy.integrate.simps(stokes[0,:] * pV, x=wavelength)		

pl.close('all')
f, ax = pl.subplots()
for j in range(4):
	ax.plot(Bpar[0,j,:], SV[0,j,:], color=pal[j], label=r'$\theta_B=${0}'.format(thetas[j]))
	ax.plot(Bpar[1,j,:], SV[1,j,:], color=pal[j])
	
ax.set_xlabel('B$_\parallel$')
ax.set_ylabel('S$_V$')
ax.set_xlim([-3500,3500])
pl.legend()

f, ax = pl.subplots()
for j in range(4):
	ax.plot(Bperp[j,:], SQ[j,:], color=pal[j], label=r'$\theta_B=${0}'.format(thetas[j]))	
	
ax.set_xlabel('B$_\perp$')
ax.set_ylabel('S$_Q$')
ax.set_xlim([0,3500])
pl.legend(loc='upper left')



    


















    
Out[2]:








<matplotlib.legend.Legend at 0x5b95bd0>










    






/usr/pkg/python/Python-2.7.6/lib/python2.7/site-packages/matplotlib-1.3.1-py2.7-linux-x86_64.egg/matplotlib/font_manager.py:1236: UserWarning: findfont: Font family ['Arial'] not found. Falling back to Bitstream Vera Sans
  (prop.get_family(), self.defaultFamily[fontext]))

Center-of-gravity





In [3]:



    


lambda0 = 5250.21
JUp = 0.0
JLow = 1.0
gUp = 0.0
gLow = 3.0
lambdaStart = 5249.5
lambdaStep = 0.015
nLambda = 100
wavelength = lambdaStart + np.arange(nLambda) * lambdaStep

lineInfo = np.asarray([lambda0, JUp, JLow, gUp, gLow, lambdaStart, lambdaStep])

s = pymilne.milne(nLambda, lineInfo)

BField = 100.0
BTheta = 0.0
BChi = 0.0
VMac = 0.0
damping = 0.04
B0 = 1.0
B1 = 12.0
mu = 1.0
VDop = 0.032
kl = 6.1

nB = 100
B = np.linspace(0.0,3500.0,nB)

Bpar = np.zeros((4,nB))
BEstimated = np.zeros((4,nB))

pal = sn.color_palette()

thetas = [15, 30, 60, 75]

for j in range(4):
	for i in range(nB):
		modelSingle = np.asarray([B[i], thetas[j], 0.0, VMac, damping, B0, B1, VDop, kl])
		stokes = s.synth(modelSingle,mu)
		
		Ic = stokes[0,0]
		
		Iplus = 0.5*(stokes[0,:] + stokes[3,:])
		Iminus = 0.5*(stokes[0,:] - stokes[3,:])
		
		lambdaPlus = scipy.integrate.simps((0.5*Ic - Iplus) * wavelength, x=wavelength) / scipy.integrate.simps((0.5*Ic - Iplus), x=wavelength)
		lambdaMinus = scipy.integrate.simps((0.5*Ic - Iminus) * wavelength, x=wavelength) / scipy.integrate.simps((0.5*Ic - Iminus), x=wavelength)
		
		BEstimated[j,i] = 1.071e9 / (3.0 * 5250.21**2) * (lambdaPlus-lambdaMinus)*1e3
		Bpar[j,i] = B[i] * np.cos(thetas[j]*np.pi/180.0)		

pl.close('all')
f, ax = pl.subplots()
for j in range(4):
	ax.plot(Bpar[j,:], (BEstimated[j,:]-Bpar[j,:]) / Bpar[j,:], color=pal[j], label=r'$\theta_B=${0}'.format(thetas[j]))	
	
ax.set_xlabel('B$_\parallel$')
ax.set_ylabel('Relative error')
ax.set_xlim([0,3500])
pl.legend()



    


















    
Out[3]:








<matplotlib.legend.Legend at 0x5b5acd0>

Influence of filling factor on inclination of magnetic field





In [15]:



    


thetaReal = np.linspace(0.0,180.0,200)
x = np.asarray([1.0,0.1,0.01])
f, ax = pl.subplots()
tanTheta = np.tan(thetaReal*np.pi/180.0)
ind = np.where(tanTheta < 0.0)
for i in range(3):
    out = np.arctan2(tanTheta, np.sqrt(x[i]))
    out[ind] += np.pi
    ax.plot(thetaReal, out*180.0/np.pi, label='x={0}'.format(x[i]))
ax.set_xlabel(r'$\theta_m$')
ax.set_ylabel(r'$\theta_\mathrm{app}$')
pl.legend(loc='upper left')



    


















    
Out[15]:








<matplotlib.legend.Legend at 0x630c790>

Line ratio technique





In [20]:



    


# First line
lambda0 = 6302.5
JUp = 1.0
JLow = 1.0
gUp = 2.5
gLow = 0.0
lambdaStart = 6301.8
lambdaStep = 0.03
nLambda = 50
wavelength = lambdaStart + np.arange(nLambda) * lambdaStep

lineInfo = np.asarray([lambda0, JUp, JLow, gUp, gLow, lambdaStart, lambdaStep])

s = pymilne.milne(nLambda, lineInfo)

VMac = 0.0
damping = 0.04
B0 = 1.0
B1 = 5.0
mu = 1.0
VDop = 0.040
kl = 6.75

pV = np.exp(-(wavelength-lambda0-60e-3)**2 / 75e-3**2)

nB = 100
B = np.linspace(0.0,3500.0,nB)

SV1 = np.zeros((4,nB))
SV2 = np.zeros((4,nB))

pal = sn.color_palette()

thetas = [0, 30, 60, 85]
for j in range(4):
    for i in range(nB):
        modelSingle = np.asarray([B[i], thetas[j], 0.0, VMac, damping, B0, B1, VDop, kl])
        stokes = s.synth(modelSingle,mu)
        SV1[j,i] = scipy.integrate.simps(stokes[3,:] * pV, x=wavelength) / scipy.integrate.simps(stokes[0,:] * pV, x=wavelength)
        
# Second line
lambda0 = 6301.5
JUp = 1.0
JLow = 1.0
gUp = 1.5
gLow = 0.0
lambdaStart = 6300.8
lambdaStep = 0.03
nLambda = 50
wavelength = lambdaStart + np.arange(nLambda) * lambdaStep

lineInfo = np.asarray([lambda0, JUp, JLow, gUp, gLow, lambdaStart, lambdaStep])

s = pymilne.milne(nLambda, lineInfo)
pV = np.exp(-(wavelength-lambda0-60e-3)**2 / 75e-3**2)
thetas = [0, 30, 60, 85]
for j in range(4):
    for i in range(nB):
        modelSingle = np.asarray([B[i], thetas[j], 0.0, VMac, damping, B0, B1, VDop, kl])
        stokes = s.synth(modelSingle,mu)
        SV2[j,i] = scipy.integrate.simps(stokes[3,:] * pV, x=wavelength) / scipy.integrate.simps(stokes[0,:] * pV, x=wavelength)
        
f, ax = pl.subplots()
for i in range(4):
    ax.plot(B, SV1[i,:] / SV2[i,:])
for i in range(4):
    ax.plot(B, SV1[i,:] / SV2[i,:], color=pal[i], label=r'$\theta_B$={0}'.format(thetas[i]))
ax.set_xlabel('B [G]')
ax.set_ylabel('Line ratio')
ax.axhline(2.5/1.5, color=pal[4])
pl.legend(loc='lower left')



    


















    
Out[20]:








<matplotlib.legend.Legend at 0x7055b50>










    
























In [ ]:

Content source: aasensio/SolarnetGranada

Similar notebooks:

notebook.community | gallery | about